12. Self-calibrated CS-pMR image reconstruction from undersampled non-Cartesian data¶
In this tutorial we will reconstruct an MRI image from radial undersampled kspace measurements. Let us denote $\Omega$ the undersampling mask, the under-sampled Fourier transform now reads $F_{\Omega}$.
We use the toy datasets available in pysap, more specifically a 2D brain slice and under-sampled Cartesian acquisition over 32 channels. We compare zero-order image reconstruction with self-calibrated multi-coil Compressed sensing reconstructions (analysis vs synthesis formulation) using the FISTA algorithm for the synthesis formulation and the Condat-Vu algorithm for the analysis formulation. The multicoil data $(y_\ell)_\ell$ is collected across multiple, say $L$, channels. The sensitivity maps $(S_\ell)_\ell$ are automically calibrated from the central portion of k-space (e.g. 5%) for all channels $\ell=1, \ldots, L$.
We remind that the synthesis formulation of the non-Cartesian CS-PMRI problem reads (minimization in the sparsifying domain):
$$ \widehat{z} = \text{arg}\,\min_{z\in C^n_\Psi} \frac{1}{2} \sum_{\ell=1}^L\|y_\ell - F_\Omega S_\ell \Psi^*z \|_2^2 + \lambda \|z\|_1 $$
and the image solution is given by $\widehat{x} = \Psi^*\widehat{z}$. For an orthonormal wavelet transform, we have $n_\Psi=n$ while for a frame we may have $n_\Psi > n$.
while the analysis formulation consists in minimizing the following cost function (min. in the image domain):
$$ \widehat{x} = \text{arg}\,\min_{x\in C^n} \frac{1}{2} \sum_{\ell=1}^L \|y_\ell - F_\Omega S_\ell x\|_2^2 + \lambda \|\Psi x\|_1 \,. $$
- Author: Chaithya G R & Philippe Ciuciu
- Date: 01/06/2021
- Target: ATSI MSc students, Paris-Saclay University
# 导入 MRI 重建运算符
# Package import
from mri.operators import NonCartesianFFT, WaveletN, WaveletUD2
'''
NonCartesianFFT : 非笛卡尔傅立叶变换(Non-Cartesian FFT),用于处理 螺旋或径向 MRI 采样 的傅立叶变换。
WaveletN : 正交小波变换(如 Symmlet-8 或 Daubechies),用于 合成公式(Synthesis Formulation) 压缩感知重建。
WaveletUD2 : 未降采样双正交小波变换(Undecimated Bi-Orthogonal Wavelet Transform),用于 分析公式(Analysis Formulation) 压缩感知重建。
'''
# 导入 MRI 处理的实用工具
from mri.operators.utils import convert_locations_to_mask, \
gridded_inverse_fourier_transform_nd
'''
convert_locations_to_mask : 将 k-space 采样坐标转换为二值掩码,用于 重建欠采样的 MRI 数据。
gridded_inverse_fourier_transform_nd : 基于网格插值的逆傅立叶变换,通常用作 MRI 重建的初始估计。
'''
# 导入 MRI 图像重建器
from mri.reconstructors import SelfCalibrationReconstructor
'''
SingleChannelReconstructor : 单通道 MRI 迭代重建器,支持 FISTA、POGM 和 Condat-Vu 等压缩感知优化算法。
'''
from mri.reconstructors.utils.extract_sensitivity_maps import get_Smaps
# 导入 MRI 样本数据
import pysap
from pysap.data import get_sample_data
# Third party import
from modopt.math.metrics import ssim
from modopt.opt.linear import Identity
from modopt.opt.proximity import SparseThreshold
'''
ssim : 结构相似性指数(SSIM),用于 评估重建图像的质量,与原始 MRI 进行比较。
Identity : 单位变换(Identity Transformation),即 不对数据做任何变换,在优化过程中常用于保持数据结构。
SparseThreshold : 稀疏阈值运算(L1 正则化 / 软阈值收缩),用于 压缩感知重建中的稀疏性约束。
'''
import numpy as np
import matplotlib.pyplot as plt
Loading input data¶
# Loading input data
cartesian_ref_image = get_sample_data('2d-pmri')
#image = pysap.Image(data=np.sqrt(np.sum(cartesian_ref_image.data**2, axis=0)))
image = np.linalg.norm(cartesian_ref_image, axis=0)
# Obtain MRI non-Cartesian radial mask
radial_mask = get_sample_data("mri-radial-samples")
kspace_loc = radial_mask.data
mask = pysap.Image(data=convert_locations_to_mask(kspace_loc, image.shape))
plt.subplot(1, 2, 1)
plt.imshow(np.abs(image), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(mask, cmap='gray')
plt.title("K-space Sampling Mask")
plt.show()
$\leadsto$ The code above is to generate the mask, the key part is
radial_mask = get_sample_data("mri-radial-samples")
where we use the package to import the mri-radial-samples already existing.
Generate the kspace¶
From the 2D brain slice and the acquisition mask, we retrospectively undersample the k-space using a non cartesian acquisition mask. We then grid the kspace to get the gridded solution as a baseline.
# Get the kspace observation values for the kspace locations
fourier_op = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
#implementation='cpu'
#implementation='finufft' # Use gpuNUFFT for speed
implementation='gpuNUFFT'
)
'''
NonCartesianFFT : 定义 非笛卡尔傅立叶变换运算符,用于处理 非规则 k-space 采样(如径向、螺旋轨迹等)。
samples=kspace_loc : 指定 k-space 采样点的坐标(即 kspace_loc)。
shape=image.shape : 设置 MRI 图像的形状(匹配输入 image 的尺寸)。
implementation='finufft' : 选择 finufft 作为计算 NUFFT(非均匀傅立叶变换) 的实现方式。
fourier_op 现在是一个 运算符对象,可以用于在非笛卡尔 k-space 进行傅立叶变换。
'''
# kspace_obs = fourier_op.op(cartesian_ref_image)
kspace_obs = fourier_op.op(cartesian_ref_image.data)
Gridded solution
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
'''
np.linspace(-0.5, 0.5, num=image.shape[0]) : 生成 均匀分布的坐标点,范围在 [-0.5, 0.5] 之间。
image.shape[0] : 设定 坐标点个数,与图像的尺寸保持一致。
np.meshgrid(grid_space, grid_space) : 生成 2D 网格坐标,用于 将非均匀采样数据重新映射到规则网格。
最终结果:grid2D 是 用于插值的 2D 坐标网格
'''
grid_soln = np.asarray([
gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
tuple(grid2D), 'linear')
for kspace_obs_ch in kspace_obs
])
'''
gridded_inverse_fourier_transform_nd() : 采用 Gridding 插值方法 进行 k-space 数据的逆傅立叶变换。
kspace_loc : MRI 非笛卡尔 k-space 采样位置(NonCartesianFFT 的 samples)。
kspace_obs : 欠采样 k-space 观测数据(fourier_op.op(image) 计算的结果)。
tuple(grid2D) : 规则网格坐标(转换为 tuple)。
'linear' : 使用 线性插值(Linear Interpolation) 进行 Gridding 重建。
最终结果:grid_soln 是 基于 Gridding 方法的 MRI 影像重建结果。
'''
image_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln)**2, axis=0)))
plt.imshow(image_rec0, cmap='gray')
base_ssim = ssim(image_rec0, image)
plt.title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 3)))
plt.show()
print('The Base SSIM is : {}'.format(base_ssim))
The Base SSIM is : 0.6173104381198481
Q2
- Observe and comment the improved image quality on the zero-filled adjoint FFT re construction when using multi-coil data as compared to single coil data.
cartesian_ref_image.data.shape
(32, 512, 512)
cartesian_ref_image.data[0].reshape(1, 512, 512).shape
(1, 512, 512)
cartesian_ref_image.data[0:2, : , :].shape
(2, 512, 512)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
'''
fourier_op_multi_coil = FFT(mask=mask, shape=image.shape,
n_coils=cartesian_ref_image.shape[0])
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image)
zero_soln_multi_coil = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim_multi_coil = ssim(zero_soln_multi_coil, image)
'''
fourier_op_multi_coil = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
#implementation='cpu'
#implementation='finufft' # Use gpuNUFFT for speed
implementation='gpuNUFFT'
)
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image.data)
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_multi_coil = np.asarray([
gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
tuple(grid2D), 'linear')
for kspace_obs_ch in kspace_obs_multi_coil
])
image_rec0_multi_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_multi_coil)**2, axis=0)))
axs[0].imshow(np.abs(image_rec0_multi_coil), cmap='gray')
base_ssim_multi_coil = ssim(image_rec0_multi_coil, image)
axs[0].set_title('Gridded solution (multi_coil) : SSIM = ' + str(np.around(base_ssim_multi_coil, 3)))
axs[0].axis('off')
'''
fourier_op_single_coil = FFT(mask=mask, shape=image.shape,
n_coils=1)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image)
zero_soln_single_coil = fourier_op_single_coil.adj_op(kspace_obs[0])
base_ssim_single_coil = ssim(zero_soln_single_coil, image)
'''
fourier_op_single_coil = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='gpuNUFFT')
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image.data[0:1, : , :].reshape(512, 512))
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_single_coil = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_single_coil,
tuple(grid2D), 'linear') # 线性插值
image_rec0_single_coil = grid_soln_single_coil
axs[1].imshow(np.abs(image_rec0_single_coil), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[1].set_title('Gridded solution (single_coil) : SSIM = ' + str(np.around(base_ssim_single_coil, 3)))
axs[1].axis('off')
axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')
plt.tight_layout()
plt.show()
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
'''
fourier_op_multi_coil = FFT(mask=mask, shape=image.shape,
n_coils=cartesian_ref_image.shape[0])
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image)
zero_soln_multi_coil = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim_multi_coil = ssim(zero_soln_multi_coil, image)
'''
fourier_op_multi_coil = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
#implementation='cpu'
#implementation='finufft' # Use gpuNUFFT for speed
implementation='gpuNUFFT'
)
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image.data)
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_multi_coil = np.asarray([
gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
tuple(grid2D), 'linear')
for kspace_obs_ch in kspace_obs_multi_coil
])
image_rec0_multi_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_multi_coil)**2, axis=0)))
axs[0].imshow(np.abs(image_rec0_multi_coil), cmap='gray')
base_ssim_multi_coil = ssim(image_rec0_multi_coil, image)
axs[0].set_title('Gridded solution (multi_coil) : SSIM = ' + str(np.around(base_ssim_multi_coil, 3)))
axs[0].axis('off')
'''
fourier_op_single_coil = FFT(mask=mask, shape=image.shape,
n_coils=1)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image)
zero_soln_single_coil = fourier_op_single_coil.adj_op(kspace_obs[0])
base_ssim_single_coil = ssim(zero_soln_single_coil, image)
'''
fourier_op_single_coil = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=2,
#implementation='cpu'
#implementation='finufft' # Use gpuNUFFT for speed
implementation='gpuNUFFT'
)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image.data[0:2, : , :])
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_single_coil = np.asarray([
gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
tuple(grid2D), 'linear')
for kspace_obs_ch in kspace_obs_single_coil
])
image_rec0_single_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_single_coil)**2, axis=0)))
axs[1].imshow(np.abs(image_rec0_single_coil), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[1].set_title('Gridded solution (double_coil) : SSIM = ' + str(np.around(base_ssim_single_coil, 3)))
axs[1].axis('off')
axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')
plt.tight_layout()
plt.show()
$\leadsto$ When using multi-coil data, the reconstrction image has a sharper contrast compared to single-coil data.
Estimate Sensitivity maps (Smaps)¶
# Obtain SMaps
Smaps, SOS = get_Smaps(
k_space=kspace_obs,
img_shape=fourier_op.shape,
samples=kspace_loc,
thresh=(0.01, 0.01), # The cutoff threshold in each kspace direction
# between 0 and kspace_max (0.5)
min_samples=kspace_loc.min(axis=0),
max_samples=kspace_loc.max(axis=0),
mode='gridding',
method='linear',
n_cpu=-1,
)
h=3;w=5;
f, axs = plt.subplots(h,w)
for i in range(h):
for j in range(w):
axs[i, j].imshow(np.abs(Smaps[3 * i + j]))
axs[i, j].axis('off')
plt.subplots_adjust(wspace=0,hspace=0)
plt.show()
! pip show pynfft2
Name: pyNFFT2 Version: 1.4.3 Summary: A pythonic wrapper around NFFT Home-page: https://github.com/ghisvail/pyNFFT.git Author: Ghislain Vaillant Author-email: ghisvail@gmail.com License: Location: /home/home/miniconda3/envs/tensor/lib/python3.11/site-packages Requires: Required-by:
Setup Fourier operators with SENSE¶
# Setup Fourier Operator with SENSE. This would initialize the
# fourier operators in the GPU.
# For this we need to specify the implementation as gpuNUFFT
# and also pass the Smaps calculated above
#fourier_implementation = 'cpu'
fourier_implementation = 'gpuNUFFT'
#fourier_implementation = 'finufft'
fourier_op_sense = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
smaps=Smaps,
implementation=fourier_implementation,
)
Q1
- Use the
gpuNUIFFTbackend instead of thefinufftone to observe a gain in computing speed when estimating the Smaps and computing the SENSE (i.e. zero-order) image.
import time
def zero_order_reconstruction (implementation_) :
time_1 = time.time()
if implementation_ == "gpuNUFFT" :
fourier_op = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
implementation='gpuNUFFT'
)
elif implementation_ == "finufft" :
fourier_op = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
implementation='finufft' # Use gpuNUFFT for speed
)
# kspace_obs = fourier_op.op(cartesian_ref_image.data)
kspace_obs = fourier_op.op(cartesian_ref_image.data)
# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln = np.asarray([
gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
tuple(grid2D), 'linear')
for kspace_obs_ch in kspace_obs
])
image_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln)**2, axis=0)))
time_2 = time.time()
time_zero_order = time_2 - time_1
base_ssim = ssim(image_rec0, image)
# Obtain SMaps
time_3 = time.time()
Smaps, SOS = get_Smaps(
k_space=kspace_obs,
img_shape=fourier_op.shape,
samples=kspace_loc,
thresh=(0.01, 0.01), # The cutoff threshold in each kspace direction
# between 0 and kspace_max (0.5)
min_samples=kspace_loc.min(axis=0),
max_samples=kspace_loc.max(axis=0),
mode='gridding',
method='linear',
n_cpu=-1,
)
time_4 = time.time()
time_SMaps = time_4 - time_3
# Setup Fourier Operator with SENSE. This would initialize the
# fourier operators in the GPU.
# For this we need to specify the implementation as gpuNUFFT
# and also pass the Smaps calculated above
#fourier_implementation = 'cpu'
time_5 = time.time()
if implementation_ == "gpuNUFFT" :
fourier_implementation = 'gpuNUFFT'
elif implementation_ == "finufft" :
fourier_implementation = 'finufft'
fourier_op_sense = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
smaps=Smaps,
implementation=fourier_implementation,
)
time_6 = time.time()
time_SENSE = time_6 - time_5
return image_rec0, time_zero_order, base_ssim, time_SMaps, Smaps, time_SENSE
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
for implementation_ in ["finufft", "gpuNUFFT"] :
# Perform reconstruction using "finufft" or "gpuNUFFT"
image_rec0, time_zero_order, base_ssim, time_SMaps, Smaps, time_SENSE = zero_order_reconstruction(implementation_)
# Define grid layout: 1 row, 2 columns (left: single image, right: multiple subplots)
fig = plt.figure(figsize=(12, 6))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2]) # Left takes 1/3, Right takes 2/3
# Left: Gridded Solution (Single Image)
ax1 = plt.subplot(gs[0])
ax1.imshow(image_rec0, cmap='gray')
ax1.set_title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 3)) + ' (' + str(implementation_) + ')')
ax1.axis('off')
ax1.text(
0.5, -0.1, # Position (centered below image)
f" time_zero_order = {time_zero_order:.3f} sec\n time_SMaps = {time_SMaps:.3f} sec\n time_SENSE = {time_SENSE:.3f} sec",
ha='center', va='top', fontsize=12, transform=ax1.transAxes
)
# Right: SMaps Subplot (Multiple Images)
h, w = 3, 5 # Define grid shape for SMaps
gs2 = gridspec.GridSpecFromSubplotSpec(h, w, subplot_spec=gs[1]) # Nested grid
for i in range(h):
for j in range(w):
ax = plt.subplot(gs2[i, j])
ax.imshow(np.abs(Smaps[3 * i + j]))
ax.axis('off')
# Adjust layout to avoid overlap
plt.tight_layout()
plt.show()
/home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5) warnings.warn( [Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers. [Parallel(n_jobs=-1)]: Done 32 out of 32 | elapsed: 0.3s finished /home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5) warnings.warn(
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done 32 out of 32 | elapsed: 0.3s finished
/home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/operators/interfaces/gpunufft.py:146: UserWarning: no pinning provided, pinning existing smaps now.
warnings.warn("no pinning provided, pinning existing smaps now.")
$\leadsto$ When using gpuNUIFFT ,
- the computation of the zero-order image is slightly faster
- the estimating the Smaps is slightly faster
- and the computation of computing SENSE is significantly slower
compared to the finufft one
FISTA optimization¶
We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost
# Setup the operators
linear_op = WaveletN(wavelet_name='sym8', nb_scale=4)
# regularizer_op = SparseThreshold(Identity(), 4e-7, thresh_type="soft")
regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")
# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
fourier_op=fourier_op_sense,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
lipschitz_cst=fourier_op_sense.impl.get_lipschitz_cst(), # add
verbose=1,
)
Lipschitz constant is 16.209672927856445
WARNING: Making input data immutable.
Synthesis formulation: FISTA optimization¶
We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = pysap.Image(data=x_final)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('FISTA Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()
- mu: 1e-08 - lipschitz constant: 16.209672927856445 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7fb2a9afdc90> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (291721,) ---------------------------------------- Starting optimization...
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:14<00:00, 6.84it/s]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 14.612145459017484 seconds ----------------------------------------
Q3
- Do the same for the CS-based (synthesis formulation) solutions (FISTA or POGM).
type(cartesian_ref_image)
pysap.base.image.Image
type(cartesian_ref_image.data)
numpy.ndarray
print(kspace_loc.shape)
(32768, 2)
from mri.reconstructors import SingleChannelReconstructor
def reconstruction (coil):
if coil == "multi_coil" :
fourier_op_sense = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
smaps=Smaps,
implementation='gpuNUFFT',
)
fourier_op = NonCartesianFFT(
samples=kspace_loc,
shape=image.shape,
n_coils=cartesian_ref_image.shape[0],
implementation='gpuNUFFT'
)
kspace_obs = fourier_op.op(cartesian_ref_image.data)
# Setup the operators
linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")
# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
fourier_op=fourier_op_sense,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
lipschitz_cst=fourier_op_sense.impl.get_lipschitz_cst(), # add
verbose=1,
)
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = pysap.Image(data=x_final)
recon_ssim = ssim(image_rec, image)
elif coil == "single_coil" :
fourier_op = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='gpuNUFFT')
kspace_obs = fourier_op.op(cartesian_ref_image.data[0])
# Setup the operators
linear_op = WaveletN(wavelet_name='sym8', nb_scale=4)
regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")
# Setup Reconstructor
reconstructor = SingleChannelReconstructor(
fourier_op=fourier_op, # the **Non Cartesian Fourier transform** (NonCartesianFFT) operator only for the sampled locations ( `kspace_loc` )
linear_op=linear_op, # Symmlet-8 wavelet transform with 4 scales for multi-scale sparse representation.
regularizer_op=regularizer_op, # Identity() < - > identity linear transformation, \lambda = 6 * 1e-7, Soft Thresholding < - > L1 norm
gradient_formulation='synthesis', # using Synthesis Formulation
verbose=1,
)
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs, # kspace_obs = fourier_op.op(image)
optimization_alg='fista', # fista
num_iterations=200,
)
recon_ssim = ssim(image_rec, image)
return image_rec, recon_ssim
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
FISTA_soln_multi_coil, FISTA_ssim_multi_coil = reconstruction("multi_coil")
axs[0].imshow(np.abs(FISTA_soln_multi_coil), cmap='gray')
axs[0].set_title('Iterative FISTA Reconstruction (multi_coil) : SSIM = ' + str(np.around(FISTA_ssim_multi_coil, 3)))
axs[0].axis('off')
FISTA_soln_single_coil, FISTA_ssim_single_coil = reconstruction("single_coil")
axs[1].imshow(np.abs(FISTA_soln_single_coil), cmap='gray')
axs[1].set_title('Iterative FISTA Reconstruction (single_coil) : SSIM = ' + str(np.around(FISTA_ssim_single_coil, 3)))
axs[1].axis('off')
axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')
plt.tight_layout()
plt.show()
WARNING: Making input data immutable. WARNING: Making input data immutable.
Lipschitz constant is 16.209671020507812 - mu: 1e-08 - lipschitz constant: 16.209671020507812 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7fb2f7cc2290> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (291721,) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [00:14<00:00, 6.97it/s]
WARNING: Making input data immutable.
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 14.347829079022631 seconds ----------------------------------------
WARNING: Making input data immutable.
Lipschitz constant is 17.830640220642092 - mu: 1e-08 - lipschitz constant: 17.830640220642092 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7fb13c585ed0> - 4 - max iterations: 200 - image variable shape: (512, 512) - alpha variable shape: (291721,) ---------------------------------------- Starting optimization...
100%|██████████| 200/200 [00:10<00:00, 19.41it/s]
- final iteration number: 200 - final log10 cost value: 6.0 - converged: False Done. Execution time: 10.306300255993847 seconds ----------------------------------------
$\leadsto$ When using multi-coil data, the quality of iterative FISTA reconstruction is better to that of single-coil data.
POGM reconstruction¶
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='pogm',
num_iterations=200,
)
image_rec = pysap.Image(data=np.abs(x_final))
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('POGM Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()
- mu: 1e-08 - lipschitz constant: 16.20966911315918 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7fb2f7a0d2d0> - 4 - max iterations: 200 - image variable shape: (1, 512, 512) ---------------------------------------- Starting optimization...
100%|██████████| 200/200 [00:28<00:00, 7.04it/s]
- final iteration number: 200 - final log10 cost value: 6.0 - converged: False Done. Execution time: 28.399563284998294 seconds ----------------------------------------
Analysis formulation: Condat-Vu reconstruction¶
#linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
linear_op = WaveletUD2(
wavelet_id=24,
nb_scale=4,
)
# regularizer_op = SparseThreshold(Identity(), 4e-7, thresh_type="soft")
regularizer_op = SparseThreshold(Identity(), 1e-9, thresh_type="soft")
# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
fourier_op=fourier_op_sense,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='analysis',
verbose=1,
)
WARNING: Making input data immutable.
Lipschitz constant is 11.979288864135743
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='condatvu',
num_iterations=200,
)
image_rec = pysap.Image(data=np.abs(x_final))
plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('Condat-Vu Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()
WARNING: <class 'mri.operators.linear.wavelet.WaveletUD2'> does not inherit an operator parent.
- mu: 1e-09 - lipschitz constant: 11.979288864135743 - tau: 0.15369599612293275 - sigma: 0.5 - rho: 1.0 - std: None - 1/tau - sigma||L||^2 >= beta/2: True - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletUD2 object at 0x7fb2f7bee610> - 4 - max iterations: 200 - number of reweights: 0 - primal variable shape: (512, 512) - dual variable shape: (2621440,) ---------------------------------------- Starting optimization...
100%|██████████| 200/200 [04:16<00:00, 1.28s/it]
- final iteration number: 200 - final cost value: 1000000.0 - converged: False Done. Execution time: 258.0937590090034 seconds ----------------------------------------